scripts/Empirical example HADS/hads_depression.R

In joakimwallmark/PolyOptimalIRT:

library(TestGardener)
library(mirt)

load(file="data/hads_dataList.rda")
load(file="data/hads_parList.rda")

# ----------- read in data -----------
titlestr  <- "hads-7-anxiety-screen"

U <- scan("data/hads.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

U <- U[, c(2, 4, 6, 8, 10, 12, 14)]
n <- 7

#  drop cases containing missing data
keepindex <- U[,7] < 4
sum(keepindex) #  N = 807
U <- U[keepindex,]

# Mirt code ---------------------------------------------------------------
colnames(U) <- paste("I", seq(1, ncol(U), 1), sep = "")

hadsDepressionGrm = mirt(data = U, model=1, itemtype = 'graded', SE=T)
#save(U, hadsDepressionGrm, file="data/hads_depression_fittedmodel_mirt.RData")
for(i in 1:ncol(U)){
    ItemPlot <- itemfit(hadsDepressionGrm,
                        group.bins=15,
                        empirical.plot = i,
                        empirical.CI = .95,
                        method = 'ML')
    print(ItemPlot)
}
itemfit(hadsDepressionGrm) # item 1 has worst fit


# TestGardener code -------------------------------------------------------
# convert U from score format to index format by adding 1
U <- U + 1

# read item/options (csv)
sheet <- read.csv("data/hads_sheet_depression.csv", header = F)

key     <- NULL # NULL for scaled items

noption  <- rep(4,n)    # number of options per item

# --------- Define the option score values for each item ---------
ind       <- c(1:n)
itemVec   <- sheet[, 1] # item strings
optionMat <- sheet[,-1] # option strings

scoreList <- list() # option scores
labelList <- list() # list of just option labels
for (item in ind){
    scoreList[[item]] <- c(0:3)
    labelList[[item]] <- optionMat[item,1:noption[item]]
}

optList <- list(itemLab=itemVec, optLab=labelList, optScr=scoreList)

scrmin  <- 0
scrmax  <- 3*n
scrrng  <- c(scrmin,scrmax)

dataList <- make.dataList(U, key, optList, scrrng)

# sum score density
ndensbasis <- 15
scoreDensity(dataList$scrvec, scrrng, ndensbasis,paste(titlestr,"- sum score"))

# percentage rank density
theta     <- dataList$percntrnk
thetaQnt  <- dataList$thetaQnt
chartList <- dataList$chartList

# ---------------  Optimal scoring: cycle of smoothing/theta estimation  ------------
#  Set number of cycles and the cell array to containing the parameter
AnalyzeResult <- Analyze(theta, thetaQnt, dataList, ncycle=10, itdisp=FALSE)

parList  <- AnalyzeResult$parList
meanHvec <- AnalyzeResult$meanHvec

#  select cycle for plotting
icycle <- 10

parListi   <- parList[[icycle]]

WfdList    <- parListi$WfdList
Qvec       <- parListi$Qvec
binctr     <- parListi$binctr
theta      <- parListi$theta
arclength  <- parListi$arclength
alfine     <- parListi$alfine

#  ----------------------------------------------------------------------------
#                   Plot surprisal curves for each test question
#  ----------------------------------------------------------------------------
#  plot both the probability and surprisal curves along with data points
Wmax <- 4
Wbinsmth.plot(binctr, titlestr, Qvec, WfdList, key, dataList, Wrng=c(0,Wmax))
Wbinsmth.plot(binctr, Qvec, WfdList, dataList, Wrng=c(0,4), saveplot=F)